Demonstration for Comparing the MSIS model versions in GEODYN

This notebook shows the output for running GEODYN with the three latest versions of the MSIS series of upper atmospheric models: MSIS_86, MSIS_00, and MSIS_2.0.

The challenges that were overcome to include these three models into versions of GEODYN are detailed [add page detailing the challenges]

Input the Parameters of the Starlette-SLR GEODYN run:

• input the parameters that make up the starlette-SLR run type in GEODYN

• further specify the runs to the different density model output.

%load_ext autoreload
%autoreload 2

##################################
# INPUT PARAMETERS:
##################################
sat = 'st'
arc = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
Accel_Status = 'acceloff'
SAT_ID = 7501001

###################################
                                  #
# CHOOSE Models to compare:       #
m1 = 'msis86'                     #
m2 = 'msis00'                     #
m3 = 'msis2'                      #
###################################



##################################
# PATH TO DENSITY MODEL RUNS:
##################################

msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/'+sat+'/results/'+msis86_model+'/'+  msis86_model+'_'+ Accel_Status +'/'

msis2_model = 'msis2'
path_to_msis2 = '/data/runs_geodyn/'+sat+'/results/'+msis2_model+'/'+  msis2_model+'_'+ Accel_Status +'/'

msis00_model = 'msis00'
path_to_msis00 = '/data/runs_geodyn/'+sat+'/results/'+msis00_model+'/'+  msis00_model+'_'+ Accel_Status +'/'

if m1 == 'msis86':
    path_to_m1 = path_to_msis86
elif m1 == 'msis00':
    path_to_m1 = path_to_msis00
elif m1 == 'msis2':
    path_to_m1 = path_to_msis2


if m2 == 'msis86':
    path_to_m2 = path_to_msis86
elif m2 == 'msis00':
    path_to_m2 = path_to_msis00
elif m2 == 'msis2':
    path_to_m2 = path_to_msis2

if m3 == 'msis86':
    path_to_m3 = path_to_msis86
elif m3 == 'msis00':
    path_to_m3 = path_to_msis00
elif m3 == 'msis2':
    path_to_m3 = path_to_msis2

import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_read_geodyn_output/')
from a_ReadGEODYN import Read_GEODYN_func

Load the simulation output for each density model

(AdjustedParams_m1,
 Trajectory_m1,
 Density_m1,
 Resids_m1)  = Read_GEODYN_func(arc,
                            sat,
                            SAT_ID,
                            grav_id,
                            local_path,
                            path_to_m1,
                            False,
                            2003,
                            'SLR')

The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/msis86/msis86_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/msis86/msis86_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/msis86/msis86_acceloff/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

(AdjustedParams_m2,
 Trajectory_m2,
 Density_m2,
 Resids_m2)  = Read_GEODYN_func(arc,
                            sat,
                            SAT_ID,
                            grav_id,
                            local_path,
                            path_to_m2,
                            False,
                            2003,
                            'SLR')

The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/msis00/msis00_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/msis00/msis00_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/msis00/msis00_acceloff/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

(AdjustedParams_m3,
 Trajectory_m3,
 Density_m3,
 Resids_m3)  = Read_GEODYN_func(arc,
                            sat,
                            SAT_ID,
                            grav_id,
                            local_path,
                            path_to_m3,
                            False,
                            2003,
                            'SLR')

The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/msis2/msis2_acceloff/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/msis2/msis2_acceloff/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/msis2/msis2_acceloff/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

Compare the Coefficient of Drag

Details: - the \(C_d\) is an adjusted parameter in GEODYN. Each iteration, it is estimated and adjusted to improve the OD simulation - \(C_d\) is treated as a constant in GEODYN (for any iteration). The factor \(C_d\) varies slightly with satellite shape and atmospheric composition. However, for any geodetically useful satellite, it may be treated as a satellite dependent constant.

• $ \bar`{A}_D = -:nbsphinx-math:frac{1}{2}` C_d \frac{A_s}{m_s} \rho_d v_r :nbsphinx-math:`bar`{r}_r $

  • \(A_D\) acceleration due to atmospheric drag force

  • \(C_d\) is the satellite drag coefficient

  • \(A_s\) is the cross-sectional area of the satellite

  • \(m_s\) is the mass of the satellite

  • \(\rho_d\)is the-density of the atmosphere at the satellite position

  • \(\bar{v}_r\) is the velocity vector of the satellite relative to the atmosphere, and

  • \(v_r\) is the magnitude of the velocity vector, vr .

Question: - What does a higher \(C_d\) actually indicate in a satellite drag scenario?

labels = list(AdjustedParams_m1[5][SAT_ID]['0CD'].keys())
val_list = []
for i in AdjustedParams_m1[5][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m1[5][SAT_ID]['0CD'][i]['CURRENT_VALUE'])



import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots


last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val

labels = list(AdjustedParams_m1[5][SAT_ID]['0CD'].keys())


fig = make_subplots(
    rows=1, cols=1,
    subplot_titles=("Time Dependent Drag Coefficient (last iter)",
                    ))
val_list = []
for i in AdjustedParams_m1[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m1[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= m1,
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )


val_list = []
for i in AdjustedParams_m2[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m2[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= m2,
                               mode='markers',
                                    marker=dict(
                                    size=8,),
                                ), row=1, col=1,
                                   )

val_list = []
for i in AdjustedParams_m3[last_iter][SAT_ID]['0CD'].keys():
    val_list.append(AdjustedParams_m3[last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= m3,
                               mode='markers',
                                    marker=dict(
                                    size=8,),
                                ), row=1, col=1,
                                   )

fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Date", row=1, col=1)


iplot(fig)

Orbit Sampled Density from Models

Details: - the density is sampled from the models along the orbit of the Starlette satellite

data_nums = 1000

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][:data_nums],
                                 y=Density_m1['rho (kg/m**3)'][:data_nums],
                                    name=m1,
                                     mode='markers',
                                    marker=dict(
                                    size=4,
                                    ),
                                   ),
                                   ],
                                   )


fig.add_trace(go.Scatter(x=Density_m2['Date'][:data_nums],
                         y=Density_m2['rho (kg/m**3)'][:data_nums],
                         name= m2,
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )

fig.add_trace(go.Scatter(x=Density_m3['Date'][:data_nums],
                         y=Density_m3['rho (kg/m**3)'][:data_nums],
                         name= m3,
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )

fig.update_layout(
    title="Density along Starlette Orbit (for 4 minutes of first day)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})

iplot(fig)

data_nums = 100
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][::data_nums],
                                 y=Density_m1['rho (kg/m**3)'][::data_nums],
                                 name = m1,
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )

fig.add_trace(go.Scatter(x=Density_m2['Date'][::data_nums],
                         y=Density_m2['rho (kg/m**3)'][::data_nums],
                         name= m2,
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )

fig.add_trace(go.Scatter(x=Density_m3['Date'][::data_nums],
                         y=Density_m3['rho (kg/m**3)'][::data_nums],
                         name= m3,
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )


fig.update_layout(
    title="Density along Starlette Orbit (Every 100th datapoint)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})

iplot(fig)

Residuals (obervation residuals)

Details: - in general, residuals are (observed minus computed) - obervation is from the observation file supplied by SLR tracking data - computed is what GEODYN is computing the location of the satellite should be -

• these observation residuals show the WHAT EVEN ARE THE RESIDUALS

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
import pandas as pd

fig = go.Figure(data=[go.Scatter(x=(Resids_m1['Date']),
                                 y=Resids_m1['Residual'].values.astype(float),
                                 name =m1 ,
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )



fig.add_trace(go.Scatter(x=(Resids_m2['Date']),
                         y=Resids_m2['Residual'].values.astype(float),
                         name= m2,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.add_trace(go.Scatter(x=(Resids_m3['Date']),
                         y=Resids_m3['Residual'].values.astype(float),
                         name= m3,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.update_layout(
    title="Observation Residuals from Final Iteration",
    yaxis_title="Residuals",
    xaxis_title="Date",
    )
fig.update_layout(legend= {'itemsizing': 'constant'})


iplot(fig)

Percent Difference in Density Output

data_nums = 3
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline


data_nums = 100

A = Density_m1['rho (kg/m**3)'][::data_nums]
B = Density_m2['rho (kg/m**3)'][::data_nums]
C = Density_m3['rho (kg/m**3)'][::data_nums]

diff1 = ((A-B)/B)*100
diff2 = ((A-C)/C)*100
diff3 = ((B-C)/C)*100


fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][::data_nums],
                                 y=diff1,
                                 name = m1+' and '+m2+' % diff',
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )

fig.add_trace(go.Scatter(x=Density_m1['Date'][::data_nums],
                                 y=diff2,
                                 name = m1+' and '+m3+' % diff',
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   )

fig.add_trace(go.Scatter(x=Density_m3['Date'][::data_nums],
                                 y=diff3,
                                 name = m2+' and '+m3+' % diff',
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   )

fig.update_layout(
    title=r"Percent Difference in MSIS models along orbit of Starlette <br> Every 100th point",
    yaxis_title='%',
    xaxis_title="Date",
    )
fig.update_layout(legend= {'itemsizing': 'constant'})

iplot(fig)

data_nums = 5000
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_m1['Date'][:data_nums],
                                 y=Density_m1['Z'][:data_nums],
                                 name = m1,
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )

fig.add_trace(go.Scatter(x=Density_m2['Date'][:data_nums],
                         y=Density_m2['Z'][:data_nums],
                         name= m2,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.add_trace(go.Scatter(x=Density_m3['Date'][:data_nums],
                         y=Density_m3['Z'][:data_nums],
                         name= m3,
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )


fig.update_layout(
#     title="Density along Starlette Orbit (Every 15th datapoint)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes( exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})

iplot(fig)

# Trajectory_m2[7501001].columns

# data_nums = 5000
# import plotly.graph_objects as go
# import numpy as np
# from plotly.offline import plot, iplot
# %matplotlib inline

# fig = go.Figure(data=[go.Scatter(x=Trajectory_m1[7501001]['Date'][:data_nums],
#                                  y=Trajectory_m1[7501001]['HEIGHT'][:data_nums],
#                                  name = m1,
#                                  mode='markers',
#                                     marker=dict(
#                                     size=3,
#                                     ),
#                                    ),
#                                    ],
#                                    )

# fig.add_trace(go.Scatter(x=Trajectory_m2[7501001]['Date'][:data_nums],
#                          y=Trajectory_m2[7501001]['HEIGHT'][:data_nums],
#                          name= m2,
#                          mode='markers',
#                          marker=dict(
#                          size=3,),
#                          ),
#                          )

# fig.add_trace(go.Scatter(x=Trajectory_m3[7501001]['Date'][:data_nums],
#                          y=Trajectory_m3[7501001]['HEIGHT'][:data_nums],
#                          name= m3,
#                          mode='markers',
#                          marker=dict(
#                          size=3,),
#                          ),
#                          )


# fig.update_layout(
# #     title="Density along Starlette Orbit (Every 15th datapoint)",
# #     yaxis_title=r'$\frac{kg}{m^3}$',
#     xaxis_title="Date",
#     )
# fig.update_yaxes( exponentformat= 'power',)
# fig.update_layout(legend= {'itemsizing': 'constant'})

# iplot(fig)